{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Laminar pipe flow - Hagen–Poiseuille solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we use icoFoam to compute a laminar pipe flow, then we compare the numerical solution with the analytical solution. <br>\n",
    "\n",
    "You will find the instructions of how to run this case in the file README.FIRST. <br>\n",
    "\n",
    "By the way, to plot the numerical solution you need to use the utility sample to get the solution at the end of the pipe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Remember, to execute a cell in ipython you need to press shift-enter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The %matplotlib magic is used to specify which matplotlib backend we want to use. \n",
    "Most of the time, in the Notebook, you will want to use the inline backend which will embed plots inside the Notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now import the libraries numpy and matplotlib. <br>\n",
    "NumPy is the fundamental package for scientific computing. <br>\n",
    "Matplotlib is used for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now use loadtxt to read the input file.  <br>\n",
    "We type np.loadtxt because loadtxt belongs to numpy which we imported with the name np.<br>\n",
    "The skiprows option is used to skip header lines.  In this case, the input file does not contain a header.<br>\n",
    "All the indormation is saved in the variable data, which is a numpy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = np.loadtxt('../postProcessing/sampleDict/108/s2_U.xy', skiprows=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now compute the velocity profile according to Hagen–Poiseuille solution. <br>\n",
    "\n",
    "$$v_{radial} = v_{max} \\left[ 1 - \\dfrac{r^2}{R_{max}^2} \\right] $$\n",
    "\n",
    "where $v_{max}$ is the maximum axial velocity obtained from the simulation, $R_{max}$ is the maximum radius, and $r$ is the local radius."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vmax = 1.4\n",
    "rmax = 0.5\n",
    "\n",
    "x = np.linspace(-1*rmax, rmax, 100)\n",
    "sol = vmax*(1 - x**2/rmax**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now plot a sample velocity profile.  <br>\n",
    "To do so, we use the plot function. <br>\n",
    "As this function belongs to matplotlib.pyplot, we append plt to plot as follows: plt.plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x78a6f235ca60>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABO3klEQVR4nO3deVyU1eIG8GcWZlhkWERWQdxxBdwQy8zC0My0TdNK81e23jZut6JbeltuVte63cWyLLPVJVNbNMxMc0NJFPeNTRAYVplhH5h5f38MTJdCZZDhzPJ8P5/5o+l95WFE5pn3nPccmSRJEoiIiIgEkYsOQERERK6NZYSIiIiEYhkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEYhkhIiIioVhGiIiISCil6ADtYTKZUFhYCG9vb8hkMtFxiIiIqB0kSUJVVRVCQ0Mhl1/8+odDlJHCwkKEh4eLjkFEREQdkJ+fj549e170/ztEGfH29gZg/mY0Go3gNERERNQeer0e4eHhlvfxi3GIMtIyNKPRaFhGiIiIHMzlplhwAisREREJxTJCREREQrGMEBERkVAsI0RERCQUywgREREJxTJCREREQrGMEBERkVAsI0RERCQUywgREREJZXUZ2blzJ6ZNm4bQ0FDIZDJs3Lix3efu2bMHSqUSMTEx1n5ZIiIiclJWl5GamhpER0dj6dKlVp1XWVmJuXPn4vrrr7f2SxIREZETs3pvmilTpmDKlClWf6GHHnoIc+bMgUKhsOpqChERETm3Ltko7+OPP0Z2djY+//xzvPrqq5c9vqGhAQ0NDZb/1uv1toxHRDZiaDJBq6tHQWUdKmoMqKpvRFV9E/T1jahpMEKC1Op4lUIOjYcbvN2V0Li7QeOhRJDGHT19PaHxUF52sy0ickw2LyNnz57Fc889h127dkGpbN+XW7x4MV566SUbJyOizlJZa8ApbRVOFulxqqgKmaXVOH+hFiVVDZCky5/fHt3USoT5eiCiuycGBXsjKkSDqGBv9OruBYWcJYXIkdm0jBiNRsyZMwcvvfQSBgwY0O7zkpOTkZSUZPlvvV6P8PBwW0QkIitJkoSs0mrsz6nArzkV+DX3Agoq6y56vFopR5ivBwK6qc1XPJqvfHiplfh9h2hoNFmunFTVN6GyzoCiynqU1xhQ3dCE08VVOF1cha0nii3neLgpEB3ugzGR/hjTuztiI3zhpe6Si75E1ElkktTxzy0ymQwbNmzAjBkz2vz/lZWV8PPzg0KhsDxnMpkgSRIUCgV+/PFHXHfddZf9Onq9Hj4+PtDpdNBoNB2NS0QdpKttxI4zJfjpZAn2ZJahosbwh2N6+nkgKliDwSHe6B/kjXB/z+YSorri4ZU6gxEFlXUoqKxDdmk1ThVV4ZRWj9PFVahvNLU6ViGXYXhPH1w3MBDXDwrCoBBvDu8QCdLe92+bfnzQaDQ4evRoq+feffdd/Pzzz1i3bh169+5tyy9PRFegSFeHTUeKsPVEMQ6cuwCj6bfPLWqlHCMi/DC6tz/GRPpjeLgPNO5uNsvioVKgX2A39AvshgkDelieN5okZJdW48C5C0jLqUBaTgUKKutwKK8Sh/Iq8dbWMwj1ccf1g4IwZWgwxvbpDjmHdIjsjtVlpLq6GpmZmZb/zsnJQUZGBvz9/REREYHk5GQUFBTg008/hVwux9ChQ1udHxgYCHd39z88T0TiVTc04YejRdiYUYC9WeWt5nsMCOqG6wcF4bqoQET39IVKKX7NRIVchv5B5isxs8dEAADOX6jFrrNl2HayGLszy1Coq8dn+87hs33nEOLjjptjQnFrbE8MDPYWnJ6IWlhdRg4cOICJEyda/rtlbse8efOwcuVKFBUVIS8vr/MSEpFNSZKEjPxKfJp6Dj8cK2o17DE60g83DgvB9VFBiOjuKTBl+/X088TsMRGYPSYCdQYj9maV4cfjxdh8rAhFunq8/0s23v8lG4NDNLhrbARuiQ2Dp4pzTIhEuqI5I12Fc0aIOl99oxHfHynCp6m5OHJeZ3m+Tw8v3BobhukxYQj3d4wC0h71jUZsP1WCDYcKsP10CRqN5l993u5K3DEyHPfE90LvAC/BKYmcS3vfv1lGiFyMrrYRK/fm4pPUXMtEVJVSjmnDQ3H32AjEhPs6/YTPCzUGrD9UgM9Sc5FbXmt5/rqoQDw6sR9G9vITmI7IebCMEFErZdUN+Gh3Dj5LPYfqhiYAQKiPO+6O74VZo8LRvZtacMKuZzJJ+OVsKT7dm4sdZ0otc2TG9e2OP03sh/i+3Z2+mBHZEssIEQEAKmoMeHd7Jj7ff84yH2RgkDcemdgXU4eFQKkQPxHVHmSXVmPZL1lYf7AATc13Do2I8MXTNwzEuH4BgtMROSaWESIXV2cwYsWeHCzbkYWq5ish0T198OjEfkgYFMRbXC/i/IVafLAzG6t/zYehyVzeJgzogeemRGFQCH//EFmDZYTIRTUZTViXfh7//OkMivXmPZ6GhGrwzOQoXNM/gMMO7VSir8fS7Zn4Yn8emkwSZDLglpgwJN0wAD39nGdiL5EtsYwQuaD0cxV4ceNxnCgyby7Z088DT98wEDdHh/JKSAflltXgHz+exqYjRQDMC779aWI/LLimD9zdFJc5m8i1sYwQuZCy6ga88cMpfJV+HgDg4+GGx6/vj7vHRkCt5BtmZzicX4nXNp/E/pwKAEBkd0/87eYhuHZgoOBkRPaLZYTIBZhMEr5Iy8M/Uk5BX2+eFzJrVDiemTzQJe+OsTVJkvDdkSK8+v0JlFSZh8AShwRh0bQhCPX1EJyOyP6wjBA5udyyGjzz9RGkNX9SHxKqwcvTh3KNjC5QVd+If/10Fh/vzYXRJKGbWom/Th2EO0eHc04O0f9gGSFyUkaThI/35GDJj6dR32iCp0qBZxIH4p74SCg4L6RLndZWIXn9ERzMqwQAjO8fgMW3DuMEV6JmLCNETiinrAZPf3UY6ecuADAvzvXGbcOdatl2R9NSDv+x5TQamkzwUimQfOMg3BUXwask5PJYRoiciCRJ+Cr9PP727XHUGozwUinw16mDMXsMhwXsRXZpNZ5ZdwQHmotiwqBAvHl7NPy9VIKTEYnDMkLkJHR1jfjrhqP4vvnW0rF9/LHkjmgOBdihlqskb6achsFoQqC3Gv+cFYOruIIruSiWESInkH6uAo+vykBBZR0UchmSJg3AQxP6cm6InTteqMPjqw4hq7QGMhnw0IS+SJo0AG5cep9cTHvfv/kvg8gOSZKED3dlY+b7+1BQWYdwfw+seygej07sxyLiAIaE+uC7x67G7DHhkCTgvR1ZuPODfSjW14uORmSXWEaI7ExNQxP+tOoQXt10EkaThOkxodj8+HjERvCWXUfiqVJi8a3D8d5dI+DtrkT6uQuY+u/d2JddLjoakd1hGSGyI1ml1Zi+dA82HSmCUi7Dy9OH4J1ZMfB2dxMdjTpoyrAQfPenqxEV7I2y6gbc9eF+fLgrGw4wQk7UZVhGiOzEj8e1mP7fPcgsqUaQRo01D47F3PhI3i3jBCIDvLDhkaswIyYURpOEVzedxGOrDqHW0CQ6GpFdYBkhEkySJLy3IwsPfJaO6oYmjOntj+8euxoje/mLjkadyEOlwD9nxeClm4dAKZfh+yNFmPl+KrQ6ziMhYhkhEsjQZMJf1h3BGymnAABz43vhi/vjEOjtLjgZ2YJMJsO8cZFY9cBY+HupcKxAj+lLd+NYgU50NCKhWEaIBKmoMeDuj/ZjXfp5KJrnh7w8fShv/3QBoyP9sfGRq9A/sBuK9Q24Y1kqUo5pRcciEoa/9YgEyC6txi3v7kFaTgW81UqsuHc05sZHio5FXSiiuye+fmQcrhnQA3WNRjz0eTre/yWLE1vJJbGMEHWxw/mVuH1ZKs6V1yLc3wPrHxmHCQN6iI5FAmjc3bBi3ijMi+8FAFj8wym88v1JmEwsJORaWEaIutCO0yW484N9qKgxYHhPH2x45Cr0D/IWHYsEUirkeGn6ULwwdRAAYMWeHDyxJgMNTUbByYi6DssIURdZf/A87v/kAOoajRjfPwCrFoxFQDe16FhkJ+4f3wfvzIqBUi7Dd4cLcd/KA6hu4K2/5BpYRoi6wIe7spG09jCaTBJmxITio3mj4aVWio5FdmZGbBhW3DsanioFdmeW4c4PUlFW3SA6FpHNsYwQ2ZAkSXjnpzN4ddNJAMD9V/fG2zNjoFLynx617ZoBPbD6gbHo3nzrL/e0IVfA34hENiJJEt5IOY13fjoLAPhL4kC8cNNgyLnRHV3G8J6++OqheIT4uCOzpBoz309FQWWd6FhENsMyQmQDkiThpe9OYNkvWQCAF28ajEcn9hOcihxJnx7dsPbBeIT7e+BceS1mLkvFufIa0bGIbIJlhKiTmUwSnt9wDCv35gIAXp0xFPdd3VtsKHJI4f6eWPtgPPoEeKGgsg4z309FZkm16FhEnY5lhKgTmUwSktcfxaq0PMhlwJI7onH32F6iY5EDC/HxwOoHx2JAkHm11js/2IesUhYSci4sI0SdRJIkvPjNMaw5kA+5DPjnrBjcPrKn6FjkBAK93bH6gXgMCtGgrLoBc5bvQ24Zh2zIebCMEHWCljkiX+zPg0wGvD0zBtNjwkTHIifi76XCF/fHYWCQN4r15kKSX1ErOhZRp2AZIbpCkiThtc0nLXNE3rxtOGbEsohQ5/P3UuHz++PQt4cXCnX1mL18H++yIafAMkJ0hZb8eBrLd+UAAF67ZRjuGBUuOBE5sx7eaqxaMBa9A7xw/kId5izfB62O65CQY2MZIboC7/+ShaXbzbfvvnTzEMyJixCciFxBoMYdXy6Is9z2O3fFflTWGkTHIuowlhGiDlr7az4W/3AKAPDclCjMGxcpNhC5lBAfD3x5/1gEadQ4U1yN+St/Ra2Be9mQY7K6jOzcuRPTpk1DaGgoZDIZNm7ceMnj169fj0mTJqFHjx7QaDSIj4/Hli1bOpqXyC6kHNPiufVHAAAPTuiDhyb0FZyIXFG4vyc+uy8OPh5uOJRXiQc/S4ehySQ6FpHVrC4jNTU1iI6OxtKlS9t1/M6dOzFp0iRs3rwZ6enpmDhxIqZNm4ZDhw5ZHZbIHuzNKsPjqw7BJAGzRoXjuclRoiORCxsQ5I2P54+Gh5sCu86WIWltBowmSXQsIqvIJEnq8E+tTCbDhg0bMGPGDKvOGzJkCGbNmoWFCxe263i9Xg8fHx/odDpoNJoOJCXqHEfP63DnB6moMRiROCQIS+eMgFLB0U4S75czpbj/k1/RaJRw99gIvDJ9KGQy7oNEYrX3/bvLf4uaTCZUVVXB39//osc0NDRAr9e3ehCJll9Ri/krf0WNwYj4Pt3xrztjWUTIbkwY0ANvz4yBTAZ8vi8P7zXvi0TkCLr8N+mSJUtQXV2NmTNnXvSYxYsXw8fHx/IID+etkiSWrrYR936chrLqBkQFe+ODuSPh7qYQHYuolWnRoVh002AAwJspp/FNRoHgRETt06Vl5Msvv8RLL72EtWvXIjAw8KLHJScnQ6fTWR75+fldmJKotYYmIx747ACySmsQ4uOOlfPHwNvdTXQsojbde1Vv3N+8MeNfvjqCfdnlghMRXV6XlZHVq1fj/vvvx9q1a5GQkHDJY9VqNTQaTasHkQgmk4S/fHUE+3Mq4K1W4uP5oxHs4y46FtElPX/jIEwZGgyD0YQHPj2AzJIq0ZGILqlLysiqVaswf/58rFq1ClOnTu2KL0nUKZb8eBrfHi6EUi7De3ePRFQwizHZP7lchn/OisGICF/o65swb8WvKKniKq1kv6wuI9XV1cjIyEBGRgYAICcnBxkZGcjLywNgHmKZO3eu5fgvv/wSc+fOxVtvvYW4uDhotVpotVrodLrO+Q6IbOSrA/l4d4d5EuDiW4fh6v4BghMRtZ+7mwIfzhuNyO6eKKisw4JP01HfaBQdi6hNVpeRAwcOIDY2FrGxsQCApKQkxMbGWm7TLSoqshQTAPjggw/Q1NSERx99FCEhIZbHE0880UnfAlHnO5Bbgec3HAUAPHZdP+43Qw7J30uFlfPHwNfTDYfzK/Hs10dwBas5ENnMFa0z0lW4zgh1pfyKWsxYugflNQZMGRqMpXNGQC7neg3kuPZmlWHuR2loMkn4S+JAPDqxn+hI5CLsdp0RIntW09CEBZ8eQHmNAUNCNXhrZjSLCDm8cX0D8LebhwAA/rHlNFKOaQUnImqNZYSomckk4ck1GTilrUJANzWWzx0FT5VSdCyiTnH32F6YF98LAJC0NgMnCrmYJNkPlhGiZm9tPY2tJ4qhUsqxfO5IhPp6iI5E1KlevGkwru4XgFqDEQs+PYCy6gbRkYgAsIwQAQBSjhVh6XbznTNv3DYMsRF+ghMRdT6lQo6lc0agd4AXCirr8NiXh9Bk5C6/JB7LCLm8zJIq/HntYQDAfVf3xi2xPQUnIrIdH083vH/PSHiqFEjNLscbKadERyJiGSHXVlXfiAc+S0eNwYixffyRPCVKdCQimxsQ5I0ld0QDAJbvysF3hwsFJyJXxzJCLstkkvDntYeRXVqDYI07/jtnBHfhJZdx47AQPDihDwDgmXVHcErLCa0kDn/zkst675cs/HiiGCqFHMvuGYmAbmrRkYi61F9uGIir+wWgrtGIBz9Lh66uUXQkclEsI+SSdp8tw5IfTwMAXp4+BDHhvmIDEQmgVMjx79mxCPP1wLnyWvx5bQZXaCUhWEbI5RTr6/HE6kOQJGDWqHDcOSZCdCQiYfy9VHj/npFQKeX46WQJlu/KFh2JXBDLCLmUJqMJj606hPIaAwaFaPDS9CGiIxEJNzTMBwtvGgwAeCPlNNLPVQhORK6GZYRcyj9/OoO0nAp4qRRYOicW7m4K0ZGI7MJdcRGYFh0Ko0nCn748hIoag+hI5EJYRshl7DhdYlnY7PXbhqNPj26CExHZD5lMhsW3DkOfAC8U6eqRtDYDJhPnj1DXYBkhl1Ckq8NTazIAAHePNX8CJKLWuqmVWHrXCKiVcuw4XYplO7NERyIXwTJCTq/JaMLjqw7hQm0jhoRq8MLUwaIjEdmtQSEavNS8w+9bP57BgVzOHyHbYxkhp/efnzPxa+4FdFMr8e5dIzhPhOgyZo0Oxy2xYTCaJDyxOoPrj5DNsYyQU0vLqcB/fj4LAHh1xlD06u4lOBGR/ZPJZHh5+hBE+HuioLIOz284yvVHyKZYRshp6Wob8eTqQzBJwK2xYZgRGyY6EpHD8HZ3w79nx0Ipl2HTkSJ8deC86EjkxFhGyClJkoTkDUdQqKtHr+6eeHnGUNGRiBxOTLgvkm4YAABY9O1xZJVWC05EzoplhJzSml/zsfmoFkq5DP++Mxbd1ErRkYgc0kPX9MW4vt1R12jE46sOoaHJKDoSOSGWEXI6WaXVeOm7EwCAP98wENHcd4aow+RyGf45KwZ+nm44XqjHki2nRUciJ8QyQk6l0WhC0poM1DUaMa5vdzx4TR/RkYgcXpDGHW/eHg0A+HB3DlKzygUnImfDMkJOZen2TBw+r4PGXYm3ZkZDLpeJjkTkFCYNDsLsMeGQJODprw5DX8/bfanzsIyQ0zicX4n//JwJAHhlxlCE+HgITkTkXF6YOthyu+9L354QHYecCMsIOYU6gxFPrc2A0SThpuEhuJnLvRN1Oi+1Em/PjIZcBnx98DxSjhWJjkROgmWEnMLrP5xEdmkNgjRqvDpjKGQyDs8Q2cKoSH88NKEvACB5/VGUVNULTkTOgGWEHN6us6X4JPUcAODN26Ph66kSnIjIuT2ZMACDQzS4UNuI577m6qx05VhGyKHp6xvxzLojAIB7xvbChAE9BCcicn4qpRz/nBUDlUKOn0+VYO2BfNGRyMGxjJBD+/v3J1HUvMpq8o1RouMQuYyBwd54OtG8Ouur359Eka5OcCJyZCwj5LB+OVOKNQfyIZMB/7g9Gp4qrrJK1JXuu7oPYiN8UdXQhOT1HK6hjmMZIYekr2/Ec1+bh2fmxUdiTG9/wYmIXI9CLsM/bh8OlVKOHadLsS6dm+lRx7CMkENavPm34ZlnJg8UHYfIZfUL9EbSJPNwzcvfn4BWx7tryHosI+Rwdp4pxao084S5N28bzuEZIsEWjO+DmHBfVNU3IXn9EQ7XkNVYRsihVNU3Inn9UQDAveMiEdenu+BERKSQy7DkDvNwzfbTpfj6YIHoSORgWEbIobyZchoFlXWI8OfwDJE96RfojacSmodrvjvOxdDIKiwj5DAO5Fbgs33mxc1ev3UYh2eI7MyC8b0xLMwH+vom7l1DVmEZIYfQ0GTEc83DMzNH9cS4fgGCExHR7ykVciy+dRgUchk2HS3C1hPFoiORg7C6jOzcuRPTpk1DaGgoZDIZNm7ceNlzduzYgREjRkCtVqNfv35YuXJlB6KSK3t3exYyS6oR0E2F528cJDoOEV3E0DAf3D++NwDgxY3HUFXfKDgROQKry0hNTQ2io6OxdOnSdh2fk5ODqVOnYuLEicjIyMCTTz6J+++/H1u2bLE6LLmms8VVeHdHJgDgbzcP4d4zRHbuyesHoFd3T2j19fjHltOi45ADkElXcA+WTCbDhg0bMGPGjIse8+yzz2LTpk04duyY5bk777wTlZWVSElJadfX0ev18PHxgU6ng0aj6WhcckAmk4Q73k9F+rkLSBgUiOVzR3FHXiIHsDezDHM+3A+ZDFj3UDxG9uLChK6ove/fNp8zkpqaioSEhFbPJSYmIjU19aLnNDQ0QK/Xt3qQa/pi/zmkn7sAL5UCL08fyiJC5CDG9QvAHSN7QpKAZ78+ioYmo+hIZMdsXka0Wi2CgoJaPRcUFAS9Xo+6urY3Vlq8eDF8fHwsj/DwcFvHJDuk1dXjjRTzJd5np0Qh1NdDcCIissZfpw5CQDcVMkuqsWxHtug4ZMfs8m6a5ORk6HQ6yyM/n9tTu6JXvj+B6oYmxIT74u64XqLjEJGVfD1VWDhtCABg6Y5M5JbVCE5E9srmZSQ4OBjFxa1v7youLoZGo4GHR9ufdNVqNTQaTasHuZZfzpRi09EiyGXA328ZCrmcwzNEjmja8BBc3S8AhiYTFn57nEvFU5tsXkbi4+Oxbdu2Vs9t3boV8fHxtv7S5KDqG41Y+I15wvP8q3pjSKiP4ERE1FEymQyvzBgKlVKOnWdKsfmoVnQkskNWl5Hq6mpkZGQgIyMDgPnW3YyMDOTl5QEwD7HMnTvXcvxDDz2E7OxsPPPMMzh16hTeffddrF27Fk899VTnfAfkdN7dnolz5bUI1rjjqebdQInIcfUO8MLDE/oCAF767jjXHqE/sLqMHDhwALGxsYiNjQUAJCUlITY2FgsXLgQAFBUVWYoJAPTu3RubNm3C1q1bER0djbfeegsffvghEhMTO+lbIGeSXVqNZb+YJ7otnDYY3dRc8p3IGTx8bV9EdvdESVUD3t56RnQcsjNXtM5IV+E6I65BkiTc/dF+7Mksx7UDe+Dje0fzVl4iJ7LzTCnmrkiDXAZ8+6erMTSMQ7DOzm7WGSFqr28PF2JPZjnUSjlevplrihA5m2sG9MBNw0NgkoC/bjwGk8nuPwtTF2EZIbtQ3dCEv286CQD408R+iOjuKTgREdnCizeZh18P51fiq3Qu20BmLCNkF/6z7SxKqhoQ2d0TD0zoIzoOEdlIkMYdTyb0BwC8kXIaulpOZiWWEbIDmSXV+Gh3DgDzpFW1UiE4ERHZ0rxxkegf2A0VNQa8vZUb6RHLCAkmSRJe+u44mkwSro8KxHVRQZc/iYgcmptCjr/dbF6Z9bN953CikPuPuTqWERJqy/Fi7DpbBpVCjoXTBouOQ0Rd5Kp+AZg6zDyZ9W9cmdXlsYyQMHUGI175/gQA4MEJfdCru5fgRETUlZ6fOggebgqk5Vbg28OFouOQQCwjJMx7v2ShoLIOYb4eeOTafqLjEFEXC/P1wJ+uM//b//umk6huaBKciERhGSEh8itqseyXLADmbcY9VJy0SuSK7h/fG72aV2b9z89nRcchQVhGSIjFP5yEocmEcX27Y8rQYNFxiEgQtVKBhTeZ54t9vDsX58prBCciEVhGqMvtzy7H5qNayGXmBZC40iqRa7suKhDj+wfAYDThtc0nRcchAVhGqEsZTRJebp60eueYCAwK4V5DRK5OJpPhxZsGQyGXYcvxYuzNKhMdiboYywh1qXXp+TheqIe3uxJ/njRAdBwishMDgrxxV1wEAODl707AyH1rXArLCHWZqvpG/GOLeevwJ67vj+7d1IITEZE9eSphADTuSpzSVmHNr9y3xpWwjFCXWbo9C2XVDegd4IW58ZGi4xCRnfHzUuHJBPMV07d+PA19PfetcRUsI9Ql8sprsaJ5/5m/3jgIKiV/9Ijoj+6J74U+PbxQXmPAf3/OFB2HugjfEahLvJ5yEgajCVf3C8D1gwJFxyEiO+WmkOPFqc23+u7J4a2+LoJlhGwu/VyF5VbeF24axFt5ieiSJkYF4poBPdBolPBmCnf1dQUsI2RTkiTh1U3mdQNmjgpHVDBv5SWiy3v+xijIZcCmo0VIP3dBdByyMZYRsqnNR7U4lFcJDzcFkngrLxG1U1SwBreP7AkA+PumE9zV18mxjJDNGJpMeCPlFADzrryBGnfBiYjIkSRNGggPNwUO5lUi5ZhWdByyIZYRspnP9p1DXkUtenirsWB8H9FxiMjBBPu4Y8E15t8dr6ecgqHJJDgR2QrLCNmErrYR/95m3oHzz5MGwEutFJyIiBzRg9f0QUA3Nc6V1+LzfedExyEbYRkhm/jv9rPQ1TViYJA37hgVLjoOETkoL7XSMt/s3z+bf6+Q82EZoU6XX1GLT/aaP8Ek3xgFhZy38hJRx80c1RP9A7uhsrYR727nQmjOiGWEOt3bW89YFjibMKCH6DhE5OCUCjmSb4wCAHy8NxeFlXWCE1FnYxmhTnWySI+NGQUAgGcnR3GBMyLqFBMHBmJMb38Ymkx456czouNQJ2MZoU71ZsopSBJw0/AQDOvpIzoOETkJmUyG56aYr46sSz+Ps8VVghNRZ2IZoU6zL7sc20+XQimX4ekbBoqOQ0ROZkSEH24YHASTBPxjC5eJdyYsI9QpJEnC6z+YFzi7c0w4IgO8BCciImf0zOSBkMuAH08Uc5l4J8IyQp1iy/FiZOSbl31//Pr+ouMQkZPqF+iNO0aalwt4I+UUl4l3EiwjdMWajCb8Y4v5qsh9V/dGoDeXfSci23lyUn+olXKk5VRgx+lS0XGoE7CM0BX7+uB5ZJXWwM/TDQ9M4LLvRGRbIT4euHdcJADz1RGTiVdHHB3LCF2R+kYj3vnJvOz7oxP7QePuJjgREbmCh6/tC427Eqe0VfjuSKHoOHSFWEboiny5Pw9FunqE+Ljj7rG9RMchIhfh66nCgxP6AgD+ufUMGo3cRM+RsYxQh9UamvDuDvPSzI9d1x/ubgrBiYjIldw7LhIB3VTILa/F1+nnRcehK8AyQh22cm8uyqoN6NXdE3eM6ik6DhG5GC+1Eg9f2w8A8O9tZ9HQZBSciDqqQ2Vk6dKliIyMhLu7O+Li4pCWlnbJ49955x0MHDgQHh4eCA8Px1NPPYX6+voOBSb7oKtrxPu/ZAMAnkzoDzcFey0Rdb274iIQ4uOOQl09vtyfJzoOdZDV7yBr1qxBUlISFi1ahIMHDyI6OhqJiYkoKSlp8/gvv/wSzz33HBYtWoSTJ0/io48+wpo1a/D8889fcXgS56Nd2dDVNaJ/YDfcHB0mOg4RuSh3NwUeu868ttHS7ZmoNTQJTkQdYXUZefvtt7FgwQLMnz8fgwcPxrJly+Dp6YkVK1a0efzevXtx1VVXYc6cOYiMjMQNN9yA2bNnX/ZqCtmv8uoGfLQ7BwCQNGkAFHJuhkdE4twxqici/D1RVm3Ayr25ouNQB1hVRgwGA9LT05GQkPDbHyCXIyEhAampqW2eM27cOKSnp1vKR3Z2NjZv3owbb7zxol+noaEBer2+1YPsx7JfslBjMGJomAaThwaLjkNELs5NIceTCearI+//kg19faPgRGQtq8pIWVkZjEYjgoKCWj0fFBQErVbb5jlz5szByy+/jKuvvhpubm7o27cvrr322ksO0yxevBg+Pj6WR3h4uDUxyYaK9fX4NPUcAODPNwyETMarIkQk3vSYMPQP7AZdXSM+3JUjOg5ZyeazDnfs2IHXXnsN7777Lg4ePIj169dj06ZNeOWVVy56TnJyMnQ6neWRn59v65jUTu9uz0RDkwmjevnh2gE9RMchIgIAKOQyJE0aAABYsTsHlbUGwYnIGkprDg4ICIBCoUBxcXGr54uLixEc3Pbl+hdffBH33HMP7r//fgDAsGHDUFNTgwceeAB//etfIZf/sQ+p1Wqo1WprolEX0OrqsSrNXAyTJg3gVREisiuJQ4IxKESDk0V6fLgrB08nDhQdidrJqisjKpUKI0eOxLZt2yzPmUwmbNu2DfHx8W2eU1tb+4fCoVCYF8fibouO5d0dmTAYTRjT2x/xfbuLjkNE1IpcLrPMHfl4Tw4u1PDqiKOwepgmKSkJy5cvxyeffIKTJ0/i4YcfRk1NDebPnw8AmDt3LpKTky3HT5s2De+99x5Wr16NnJwcbN26FS+++CKmTZtmKSVk/wor67C6+arIkwn9eVWEiOzSDYODMDhEgxqDEct3ZYuOQ+1k1TANAMyaNQulpaVYuHAhtFotYmJikJKSYpnUmpeX1+pKyAsvvACZTIYXXngBBQUF6NGjB6ZNm4a///3vnfddkM21XBWJ6+2PcX0DRMchImqTTGa+OvLAZ+n4ZG8u7h/fB/5eKtGx6DJkkgOMlej1evj4+ECn00Gj0YiO43IKKutw7T+2o9EoYfUDYzG2D4doiMh+SZKEaf/djWMFejx8bV88OzlKdCSX1d73b67hTZe1dHsmGo0S4vt0ZxEhIrsnk8nw5PXmO2s+2ZuL8uoGwYnoclhG6JLOX6jFVwd+mytCROQIrh8UiGFhPqg1GPEB547YPZYRuqSl27PQaJQwrm93xPGqCBE5iJa5IwDw6d5zKOPVEbvGMkIXVVBZh3Xp5qsiTzUvJkRE5CiuiwpEdE8f1DUauSqrnWMZoYtatiPLMldkdKS/6DhERFaRyWSWHX0/S83luiN2jGWE2lSsr8ea5rkij13fT3AaIqKOuX5QoGXdkY/38OqIvWIZoTa9/0s2DM170MRzrggROSjz1RHzB6qP9+RCV8cdfe0Rywj9QVl1A75MM+/M+9j1XG2ViBxb4pBgDAjqhqqGJnyyN1d0HGoDywj9wfJd2ahvNCG6pw+u6c/VVonIscnlMjw60Xx1ZMWeHFQ3NAlORL/HMkKtXKgx4LPU5qsi1/GqCBE5h5uGh6JPgBcqaxstv+PIfrCMUCsr9uSg1mDE4BANrh8UKDoOEVGnUMhleKT56siHu7JRa+DVEXvCMkIWurpGrNyTCwB47Lp+vCpCRE5lekwowv09UF5jwJf780THof/BMkIWn+7NRVVDEwYEdUPikGDRcYiIOpWbQo5HrjVfHflgZzbqG42CE1ELlhECANQamvBx8yzzRyf2g1zOqyJE5HxuG9ETIT7uKKlqwPqDBaLjUDOWEQIArE7LR0WNARH+npg6LER0HCIim1Ap5Vgwvg8AYNkvWWgymgQnIoBlhAAYmkxY3ryr5YMT+kCp4I8FETmvO8eEw8/TDXkVtdh0tEh0HALLCAHYeKgARbp6BHqrcduInqLjEBHZlKdKif+7qjcA4L0dWZAkSXAiYhlxcUaThGW/ZAEA7h/fG+5uCsGJiIhsb258JLxUCpzSVuHnUyWi47g8lhEXl3JMi+yyGvh4uGFOXC/RcYiIuoSPpxvujjf/zlu6PZNXRwRjGXFhkiTh3R2ZAIB54yLRTa0UnIiIqOvcd3VvqJRyHMyrxP6cCtFxXBrLiAv75Uwpjhfq4alSYP64SNFxiIi6VKC3O2aOMs+TW7o9U3Aa18Yy4sLe3WGeKzJ7TAT8vFSC0xARdb0Hr+kLhVyGXWfLcPS8TnQcl8Uy4qLSz11AWk4F3BQy3D++t+g4RERChPt74uboUADAsp1ZgtO4LpYRF/VB8z+6W2LDEOLjITgNEZE4D04wL4L2w9EinCuvEZzGNbGMuKCs0mr8eKIYAPDANX0EpyEiEisqWINrB/aASQI+3JUjOo5LYhlxQR/uyoYkAQmDgtAv0Ft0HCIi4R68pi8AYO2BfJRXNwhO43pYRlxMSVU9vk43bw710AReFSEiAoCxffwR3dMHDU0mfJJ6TnQcl8My4mJW7smFwWjCyF5+GBXpLzoOEZFdkMlkeHCC+erIp6m5qDU0CU7kWlhGXEh1QxM+22du/A9yrggRUSuJQ4LRq7snKmsbsfbXfNFxXArLiAtZtT8PVfVN6NvDCwmDgkTHISKyKwq5DAvGmz+oLd+Vg0ajSXAi18Ey4iIMTSZ8tNs8S/zBa/pCLpcJTkREZH9uH9kTAd1UKKisw+ajRaLjuAyWERfx7eFCaPX1CPRWY3psqOg4RER2yd1NgXubt8dY9ks2N9DrIiwjLkCSJHy4KxsAMP+q3lArFYITERHZr7vH9oKnSoGTRXrsySwXHcclsIy4gF1ny3BKWwVPlQJz4iJExyEismu+nirMHBUOAFje/EGObItlxAV82DxXZOaocPh4uAlOQ0Rk/+ZfFQmZzLy7+ZniKtFxnB7LiJM7ra3CzjOlkMuA+67mhnhERO3Rq7sXEgcHAwA+4hLxNscy4uRa5opMHhqMcH9PwWmIiBzHgmvMH+A2HCpAaRWXiLelDpWRpUuXIjIyEu7u7oiLi0NaWtolj6+srMSjjz6KkJAQqNVqDBgwAJs3b+5QYGq/kqp6fJNRCAC4fzwXOSMissbIXv6IjfCFwWjCZ6m5ouM4NavLyJo1a5CUlIRFixbh4MGDiI6ORmJiIkpKSto83mAwYNKkScjNzcW6detw+vRpLF++HGFhYVccni7t073nLEu/j4jwEx2HiMjhtCyC9tm+c6gzGAWncV5Wl5G3334bCxYswPz58zF48GAsW7YMnp6eWLFiRZvHr1ixAhUVFdi4cSOuuuoqREZGYsKECYiOjr7i8HRxtYYmfL7fvPT7gvGcK0JE1BGJQ4IR7u+BC7WN+PrgedFxnJZVZcRgMCA9PR0JCQm//QFyORISEpCamtrmOd9++y3i4+Px6KOPIigoCEOHDsVrr70Go/HiDbOhoQF6vb7Vg6zzdfp5VNY2old3T0xqnoRFRETWUchl+L+rzB/oVuzOgcnERdBswaoyUlZWBqPRiKCg1vuaBAUFQavVtnlOdnY21q1bB6PRiM2bN+PFF1/EW2+9hVdfffWiX2fx4sXw8fGxPMLDw62J6fJMJsmy9Pv/XdUbCi79TkTUYTNHhcPbXYnsshpsO9X2lAS6Mja/m8ZkMiEwMBAffPABRo4ciVmzZuGvf/0rli1bdtFzkpOTodPpLI/8fO6eaI2fT5Ugt7wWGnclbh/ZU3QcIiKH5qVWWhaMXLGbt/naglVlJCAgAAqFAsXFxa2eLy4uRnBw20MBISEhGDBgABSK35YgHzRoELRaLQwGQ5vnqNVqaDSaVg9qvxV7zP9YZsdFwEutFJyGiMjxzYuPhEIuQ2p2OU4UcupAZ7OqjKhUKowcORLbtm2zPGcymbBt2zbEx8e3ec5VV12FzMxMmEy/bcV85swZhISEQKVSdTA2XczJIj32ZpVDIZdhbnyk6DhERE4h1NcDU4aaP3R/vIdXRzqb1cM0SUlJWL58OT755BOcPHkSDz/8MGpqajB//nwAwNy5c5GcnGw5/uGHH0ZFRQWeeOIJnDlzBps2bcJrr72GRx99tPO+C7Jo+UcyeWgwwnw9BKchInIe/9e8ivU3GYUoq+YiaJ3J6mv4s2bNQmlpKRYuXAitVouYmBikpKRYJrXm5eVBLv+t44SHh2PLli146qmnMHz4cISFheGJJ57As88+23nfBQEAyqsbsLF5kbOW2d9ERNQ5RkT4ISbcFxn5lfhyfx4ev76/6EhOQyZJkt3fp6TX6+Hj4wOdTsf5I5fwn21n8dbWM4gO98XGR8ZBJuNdNEREnenbw4V4fNUh9PBWY/ezE6FWKi5/kgtr7/s396ZxEoYmEz7dZ17k7P+uimQRISKygSlDgxGscUdpVQM2HSkSHcdpsIw4iU1HC1Fa1YAgjRo3DgsRHYeIyCm5KeSYO64XAOCj3TlwgMEFh8Ay4gQk6bdFzubGR8JNwb9WIiJbmT06Au5uchwv1OPX3Aui4zgFvms5gQPnLuBYgR5qpRyzx0SIjkNE5NT8vFS4Jda8oCQXQescLCNOoOV23ltHhMHfi2u3EBHZ2v9dFQkA+PGEFvkVtWLDOAGWEQdXWFmHLcfNK+LeO4638xIRdYX+Qd4Y3z8AJgmWHdKp41hGHNwX+8/BaJIQ36c7BgZ7i45DROQy5jWvcr3m13zUGS6+Ez1dHsuIA6tvNGJVmnkTwXnjIsWGISJyMROjAhHu74HK2kZ8k1EgOo5DYxlxYN8dLkRFjQFhvh5IGBQoOg4RkUtRyGWYOzYSALByby5v870CLCMOSpIkfJKaCwC4e2wvKHk7LxFRl5s5Khwebgqc0lYhLadCdByHxXcwB3Uw77fbee8cHS46DhGRS/LxdMOM2DAAsHxAJOuxjDiolXvNs7enx4TCj7fzEhEJM695RdYtx4tRWFknOI1jYhlxQMX6evxw1LwnAieuEhGJFRWswdg+/jCaJHzB23w7hGXEAX2xPw9NJgmjI/0wJNRHdBwiIpd3b/MHw1Vp+ahv5G2+1mIZcTCGJhO+3J8HgFdFiIjsRcKgIIT6uKOixoDvuZuv1VhGHMwPx4pQVm3enTdxSLDoOEREBECpkOPuePPckU85kdVqLCMO5vN95vHIOWN6cXdeIiI7MmtUOFQKOY6c1+FwfqXoOA6F72YO5GSRebtqpVyG2WN4Oy8RkT3p3k2NqcNDAACf7eNEVmuwjDiQlh/uxKHBCNS4C05DRES/d0/zUM13hwtxocYgOI3jYBlxEPr6Rmw8ZN774J6xvQSnISKitsSG+2JIqAYNTSasSz8vOo7DYBlxEOvTz6PWYET/wG6I6+0vOg4REbVBJpNZPjB+vv8cTCbuV9MeLCMOQJIkyxDNPfG9IJPJBCciIqKLuTkmFN7uSpwrr8XOs6Wi4zgElhEHkJpdjqzSGnipFLileQ8EIiKyT54qJe4Yab7J4HNOZG0XlhEH0PLDfMuIMHi7uwlOQ0REl3PX2AgAwLZTJcivqBWcxv6xjNg5ra4eW44XAwDu5sRVIiKH0LdHN1zdLwCSBHyZlic6jt1jGbFzq9LyYDRJGBPpj6hgjeg4RETUTi0fINf8mo+GJu5XcyksI3as0WjCquZG3bLMMBEROYaEQYEIad6vJuWYVnQcu8YyYse2nSxBSVUDunupMJn70BARORSlQo47R5vnjnyxj0M1l8IyYsdaxhlnjg6HSsm/KiIiRzNrdDgUchnScitwtrhKdBy7xXc4O5VXXoudZ8z3p89ubtZERORYgn3ccX1UIADgi/28OnIxLCN2quWqyDUDeiCiu6fgNERE1FF3NU9k/frgedQZOJG1LSwjdsjQZMJXB/IBAHfF8aoIEZEjG98vAOH+Hqiqb8L3RwpFx7FLLCN2aMtxLcprDAjSqC2X94iIyDHJ5TLMGWO+OsKhmraxjNihL/abV1ydNToCSgX/ioiIHN0do3rCTSFDRn4ljhfqRMexO3ynszOZJdXYl10BuQy4c3S46DhERNQJArqpkdi8RMOXvDryBywjdqZlkbProoIQ6ushOA0REXWWu+LMQzUbDxWguqFJcBr7wjJiR+objViXfh4AJ64SETmbsX380aeHF2oMRnyTUSA6jl3pUBlZunQpIiMj4e7ujri4OKSlpbXrvNWrV0Mmk2HGjBkd+bJOb/PRIujqGhHm64FrBvQQHYeIiDqRTCbDnDHmD5ocqmnN6jKyZs0aJCUlYdGiRTh48CCio6ORmJiIkpKSS56Xm5uLp59+GuPHj+9wWGe3Os18O++dzSv2ERGRc7ltRE+olHIcL9Tj6HlOZG1hdRl5++23sWDBAsyfPx+DBw/GsmXL4OnpiRUrVlz0HKPRiLvuugsvvfQS+vTpc0WBnVVmSRXSciugkMtwxyhOXCUickZ+XipMGdo8kTWNV0daWFVGDAYD0tPTkZCQ8NsfIJcjISEBqampFz3v5ZdfRmBgIO677752fZ2Ghgbo9fpWD2fXclVk4sBABPu4C05DRES20rJ53rcZBajhRFYAVpaRsrIyGI1GBAUFtXo+KCgIWm3b2yPv3r0bH330EZYvX97ur7N48WL4+PhYHuHhzn2loKHJiK8Pmieuzh7j3N8rEZGrG9vHH70DzBNZuSKrmU3vpqmqqsI999yD5cuXIyAgoN3nJScnQ6fTWR75+fk2TCnej8eLcaG2EcEad0zgxFUiIqcmk8kwq3kdqVVpzv3+1l5Kaw4OCAiAQqFAcXFxq+eLi4sRHBz8h+OzsrKQm5uLadOmWZ4zmUzmL6xU4vTp0+jbt+8fzlOr1VCr1dZEc2gta4vMHNWTK64SEbmA20b0xJItp5GRX4mTRXoMCtGIjiSUVe98KpUKI0eOxLZt2yzPmUwmbNu2DfHx8X84PioqCkePHkVGRoblcfPNN2PixInIyMhw+uGX9jhXXoO9WeWQyYCZXHGViMgl9PBW44Yh5ikPqzmR1borIwCQlJSEefPmYdSoURgzZgzeeecd1NTUYP78+QCAuXPnIiwsDIsXL4a7uzuGDh3a6nxfX18A+MPzrmr1r+ZLdNf074Gefp6C0xARUVe5c3QENh/VYsOhAiTfOAjubgrRkYSxuozMmjULpaWlWLhwIbRaLWJiYpCSkmKZ1JqXlwe5nEMN7dFoNOGrA5y4SkTkiq7uF4Cefh44f6EOm48W4dYRPUVHEkYmSZIkOsTl6PV6+Pj4QKfTQaNxnnG1lGNFeOjzgwjopkZq8nVw43wRIiKX8p9tZ/HW1jMYHemHrx4aJzpOp2vv+zff/QRqmUV9+8ieLCJERC7ojlHhkMuAX3MvILOkSnQcYfgOKEhhZR12ni0FYF7+nYiIXE+wjzuuiwoEAKxtHrZ3RSwjgqxLPw9JAuJ6+yMywEt0HCIiEmRm8xYg6w+eR6PRJDiNGCwjAphMEtYeMA/RzOJVESIilzYxKhAB3dQoqzZg28lLbzrrrFhGBEjNLsf5C3XwVisxZWiI6DhERCSQm0KO20aGAYDlg6qrYRkRYE3z2iI3x4TCQ+W695UTEZFZy1DNjtMlKNbXC07T9VhGupiuthEpx82bCnKIhoiIAKBvj24YHekHk2SeU+hqWEa62MaMAhiaTIgK9sawMB/RcYiIyE60XB356kA+HGAJsE7FMtLFWsYDZ44Kh0wmE5yGiIjsxY3DQuClUiC3vBb7cypEx+lSLCNd6FiBDscL9VAp5LglNkx0HCIisiNeaiWmRYcCcL2JrCwjXajlh2vSkCD4eakEpyEiInvTsnv75qNF0Nc3Ck7TdVhGukh9oxEbDxUAAGaN4sRVIiL6o9hwX/QP7Ib6RhO+O1woOk6XYRnpIluOa6Gvb0KYrweu6hcgOg4REdkhmUxmudNy7a+uM1TDMtJFWm7Vum1kTyjknLhKRERtmxEbBqVchsPndThT7Bqb57GMdIGCyjrsziwDANw+oqfgNEREZM8Cuqktm+e5ypojLCNdYP3/bIoX0d1TdBwiIrJzt480f3Bdf7DAJTbPYxmxMUmSsO6gudnewYmrRETUDubN81Qoq27AL6dLRcexOZYRG/s19wLOldfCS6XAjcOCRcchIiIH4KaQY0aMeT0qVxiqYRmxsXXp5tnQNw4LgadKKTgNERE5ittHmYdqtp0qRkWNQXAa22IZsaFaQxM2HSkCwCEaIiKyTlSwBsPCfNBolPBNRoHoODbFMmJDPxzVosZgRK/unhgd6Sc6DhEROZiWiaxfHXDuoRqWERv6qnmI5vYRPbkpHhERWe3m6FCoFHKcKNLjeKFOdBybYRmxkbzyWuzLroBMZl7ojIiIyFp+XiokDHb+NUdYRmzk6+bbea/uF4BQXw/BaYiIyFHdMdI85/CbjEIYmpxzzRGWERswmSRLGbmdV0WIiOgKjO8fgEBvNSpqDPj5VLHoODbBMmIDabkVOH+hDt5qJRKHcG0RIiLqOKVCjltGmNcc+fqgc95VwzJiA183j+tNHR4CdzeF4DREROTobmve12z7qRKUVzcITtP5WEY6WZ3BiM1HzWuLcOIqERF1hgFB3hgW5oMmk4TvDheKjtPpWEY62Zbj5rVFIvw9MaoX1xYhIqLOcZsTD9WwjHSylomrt44I49oiRETUaaZFh0Ipl+FogQ5niqtEx+lULCOdSKurx57MMgDArbEcoiEios7TvZsaE6PMa460fPB1FiwjnWhjRgFMEjAm0h8R3T1FxyEiIifTMlSz8VABjCZJcJrOwzLSSSRJstxFc2vzDwsREVFnmhgVCF9PNxTrGyxX4p0By0gnOVagx9mSaqiVctw4PER0HCIickJqpQLThocCANY70VANy0gnaRm/u2FIMDTuboLTEBGRs2pZNiLluBZV9Y2C03QOlpFOYGgy4dvm+75v4xANERHZUHRPH/Tp4YX6RhN+OKYVHadTsIx0gh2nS1BRY0APbzWu7hcgOg4RETkxmUxmWZH1ayfZybdDZWTp0qWIjIyEu7s74uLikJaWdtFjly9fjvHjx8PPzw9+fn5ISEi45PGOaMMh8wI0M2JCoVSw3xERkW3dEhsGmQzYn1OB8xdqRce5Yla/c65ZswZJSUlYtGgRDh48iOjoaCQmJqKkpKTN43fs2IHZs2dj+/btSE1NRXh4OG644QYUFDjHCnK62kZsO2n+3m/h2iJERNQFQn09MLZ3dwDANxmOvzy81WXk7bffxoIFCzB//nwMHjwYy5Ytg6enJ1asWNHm8V988QUeeeQRxMTEICoqCh9++CFMJhO2bdt2xeHtweZjRTAYTYgK9sbgUI3oOERE5CJadvJdf/A8JMmx1xyxqowYDAakp6cjISHhtz9ALkdCQgJSU1Pb9WfU1taisbER/v7+Fz2moaEBer2+1cNebWjeI2BGLCeuEhFR15kyNBhqpRxZpTU4Xmi/75PtYVUZKSsrg9FoRFBQUKvng4KCoNW2b0bvs88+i9DQ0FaF5vcWL14MHx8fyyM8PNyamF0mv6IWabkVkMmA6TGhouMQEZEL8XZ3w6TB5vfj9Q6+eV6XzrZ8/fXXsXr1amzYsAHu7u4XPS45ORk6nc7yyM/P78KU7fdNhvkvP75Pd4T4eAhOQ0REruaW5qvy3x4uRJPRJDhNxymtOTggIAAKhQLFxcWtni8uLkZwcPAlz12yZAlef/11/PTTTxg+fPglj1Wr1VCr1dZE63KSJP12Fw2HaIiISIBrBvSAv5cKZdUN2JNVjgkDeoiO1CFWXRlRqVQYOXJkq8mnLZNR4+PjL3rem2++iVdeeQUpKSkYNWpUx9PakWMFemSV1kCtlGPK0EsXMSIiIltwU8gxrXkLkg0OvDy81cM0SUlJWL58OT755BOcPHkSDz/8MGpqajB//nwAwNy5c5GcnGw5/o033sCLL76IFStWIDIyElqtFlqtFtXV1Z33XQiw/pD5L33S4CB4c/l3IiISpOXq/JbjxahpaBKcpmOsGqYBgFmzZqG0tBQLFy6EVqtFTEwMUlJSLJNa8/LyIJf/1nHee+89GAwG3H777a3+nEWLFuFvf/vblaUXpMlownfNy7/fwiEaIiISKCbcF70DvJBTVoMtx7W4dYTjrXklkxzg5mS9Xg8fHx/odDpoNOLX8th+ugTzP/4V/l4q7H/+erhx1VUiIhLoXz+dxT9/OoPx/QPw2X1xouNYtPf9m++iHbCxeeLqtOEhLCJERCTcjFjz8hJ7MstQoq8XnMZ6fCe1UnVDE7YcN6+pwrtoiIjIHvTq7oUREb4wSbDsIu9IWEastPWEFvWNJkR290RMuK/oOERERAB+m8O4McPxFkBjGbHSxkPmxjk9JgwymUxwGiIiIrOpw0OhlMtwrECPzBLHumOVZcQKZdUN2J1ZBoDLvxMRkX3x91JhfP8AAMC3DnZ1hGXECpuOFMFokjC8pw/69OgmOg4REVErLXMZvzlc6FA7+bKMWKFlL5qbo3lVhIiI7E/CoCB4uClwrrwWh8/rRMdpN5aRdsorr8XBvErIZCwjRERkn7zUStwwxLwIacsyFI6AZaSdvj1s/ksd17c7AjUX33GYiIhIpJY5jd8fKXKYnXxZRtpBkiRszPjtLhoiIiJ7Nb7/bzv57s0qFx2nXVhG2uFEkfk2KZVSjsncoZeIiOyYm0KOqcPMO/l+k+EYC6CxjLTDt81/mddHBULDHXqJiMjOtQzVbDmuRX2jUXCay2MZuQyTSbIsrcshGiIicgQje/mhp58HqhuasO1kieg4l8UychlpuRUo0tXD212Jawf2EB2HiIjosmQymeXOT0dYHp5l5DJaxttuHBoCdzeF4DRERETt07IA2o7TJdDVNgpOc2ksI5dgaDLhh2NFAICbufw7ERE5kAFB3ogK9kajUULK8SLRcS6JZeQSdmeWorK2ET281Rjbp7voOERERFaZ1jxU0zL30V6xjFzCd4fNTXLqsBAo5Nyhl4iIHEvLvJHUrHKUVNULTnNxLCMXUWcw4sfjWgC/NUsiIiJHEu7viZhwX5gkYPMR+x2qYRm5iO2nS1BjMCLM1wMjInxFxyEiIuqQlqsj37GMOJ6Whc6mRYdCJuMQDREROaapw0MgkwHp5y7g/IVa0XHaxDLSBn19I34+bV4khjv0EhGRIwvSuGNsb/NNGC1zIe0Ny0gbth4vhqHJhL49vDAoxFt0HCIioivSMvfxOzu9q4ZlpA0tt0DdHB3GIRoiInJ4U4YGQymXWTZ+tTcsI79TUWPA7swyAMC06BDBaYiIiK6cn5cK4/sHALDPqyMsI7+z+WgRjCYJQ8M06NOjm+g4REREnaJlJfHvDhdCkiTBaVpjGfmdliGaacM5cZWIiJzHpMHBUCvlyC6rwfFCveg4rbCM/A+trh6/5lYAAG7iXTREROREuqmVuH5QIAD7Wx6eZeR/bDpaBEkCRvXyQ5ivh+g4REREnarlqv+mI0V2NVTDMvI/vj9iboo3DefEVSIicj4TowLhqVKgoLIOh/IrRcexYBlpdv5CLQ7lVUImA24cxjJCRETOx91NgYRBQQCA7+1oATSWkWabmtfsHxPpj0CNu+A0REREttFy9X/z0SKYTPYxVMMy0mzTUXMZ4cRVIiJyZhMG9oC3Wgmtvh7peRdExwHAMgIAOFdegyPndZDLzKvUEREROSu1UoFJQ1qGauzjrhqWEQDfNw/RjOsbgIBuasFpiIiIbMsyVHNMC6MdDNWwjOC3MsK7aIiIyBVc3a8HfDzcUFrVgLScCtFxWEaySqtxskgPpVyGxCEcoiEiIuenUsqR2DJUc0T8UE2HysjSpUsRGRkJd3d3xMXFIS0t7ZLHf/XVV4iKioK7uzuGDRuGzZs3dyisLbTcRXNVvwD4eakEpyEiIuoaNzUvgJZyTIsmo0loFqvLyJo1a5CUlIRFixbh4MGDiI6ORmJiIkpKSto8fu/evZg9ezbuu+8+HDp0CDNmzMCMGTNw7NixKw7fGbjQGRERuaL4vt3h5+mG8hoD9mWLHaqxuoy8/fbbWLBgAebPn4/Bgwdj2bJl8PT0xIoVK9o8/l//+hcmT56Mv/zlLxg0aBBeeeUVjBgxAv/973+vOPyVOlNchTPF1VAp5LiBQzRERORC3BRyTB5q/iAueqjGqjJiMBiQnp6OhISE3/4AuRwJCQlITU1t85zU1NRWxwNAYmLiRY8HgIaGBuj1+lYPW2i5pemaAQHw8XCzydcgIiKyV9OaRwVSjmvRKHCoxqoyUlZWBqPRiKCgoFbPBwUFQavVtnmOVqu16ngAWLx4MXx8fCyP8PBwa2K2iyRJ+L5lobPhXOiMiIhcT1yf7gjopkJlbSP2ZJYJy2GXd9MkJydDp9NZHvn5+Z3+NSQJeCYxCtNjQi1bKhMREbkShVyGB67pg2cnRyEqWCMsh9KagwMCAqBQKFBcXNzq+eLiYgQHtz3nIjg42KrjAUCtVkOttu3iY3K5DJOHBmMyV1wlIiIX9sA1fUVHsO7KiEqlwsiRI7Ft2zbLcyaTCdu2bUN8fHyb58THx7c6HgC2bt160eOJiIjItVh1ZQQAkpKSMG/ePIwaNQpjxozBO++8g5qaGsyfPx8AMHfuXISFhWHx4sUAgCeeeAITJkzAW2+9halTp2L16tU4cOAAPvjgg879ToiIiMghWV1GZs2ahdLSUixcuBBarRYxMTFISUmxTFLNy8uDXP7bBZdx48bhyy+/xAsvvIDnn38e/fv3x8aNGzF06NDO+y6IiIjIYckkSRK/Q85l6PV6+Pj4QKfTQaMRN8GGiIiI2q+97992eTcNERERuQ6WESIiIhKKZYSIiIiEYhkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEYhkhIiIioVhGiIiISCirl4MXoWWRWL1eLzgJERERtVfL+/blFnt3iDJSVVUFAAgPDxechIiIiKxVVVUFHx+fi/5/h9ibxmQyobCwEN7e3pDJZKLjCKXX6xEeHo78/Hzu02NjfK27Bl/nrsHXuWvwdW5NkiRUVVUhNDS01Sa6v+cQV0bkcjl69uwpOoZd0Wg0/EHvInytuwZf567B17lr8HX+zaWuiLTgBFYiIiISimWEiIiIhGIZcTBqtRqLFi2CWq0WHcXp8bXuGnyduwZf567B17ljHGICKxERETkvXhkhIiIioVhGiIiISCiWESIiIhKKZYSIiIiEYhlxABUVFbjrrrug0Wjg6+uL++67D9XV1e06V5IkTJkyBTKZDBs3brRtUAdn7etcUVGBxx57DAMHDoSHhwciIiLw+OOPQ6fTdWFqx7B06VJERkbC3d0dcXFxSEtLu+TxX331FaKiouDu7o5hw4Zh8+bNXZTUsVnzOi9fvhzjx4+Hn58f/Pz8kJCQcNm/FzKz9ue5xerVqyGTyTBjxgzbBnRALCMO4K677sLx48exdetWfP/999i5cyceeOCBdp37zjvvuPwS+u1l7etcWFiIwsJCLFmyBMeOHcPKlSuRkpKC++67rwtT2781a9YgKSkJixYtwsGDBxEdHY3ExESUlJS0efzevXsxe/Zs3HfffTh06BBmzJiBGTNm4NixY12c3LFY+zrv2LEDs2fPxvbt25Gamorw8HDccMMNKCgo6OLkjsXa17lFbm4unn76aYwfP76LkjoYiezaiRMnJADSr7/+annuhx9+kGQymVRQUHDJcw8dOiSFhYVJRUVFEgBpw4YNNk7ruK7kdf5fa9eulVQqldTY2GiLmA5pzJgx0qOPPmr5b6PRKIWGhkqLFy9u8/iZM2dKU6dObfVcXFyc9OCDD9o0p6Oz9nX+vaamJsnb21v65JNPbBXRKXTkdW5qapLGjRsnffjhh9K8efOk6dOnd0FSx8IrI3YuNTUVvr6+GDVqlOW5hIQEyOVy7N+//6Ln1dbWYs6cOVi6dCmCg4O7IqpD6+jr/Hs6nQ4ajQZKpUNs+2RzBoMB6enpSEhIsDwnl8uRkJCA1NTUNs9JTU1tdTwAJCYmXvR46tjr/Hu1tbVobGyEv7+/rWI6vI6+zi+//DICAwN51fQS+BvTzmm1WgQGBrZ6TqlUwt/fH1qt9qLnPfXUUxg3bhymT59u64hOoaOv8/8qKyvDK6+80u4hNFdQVlYGo9GIoKCgVs8HBQXh1KlTbZ6j1WrbPL69fw+uqCOv8+89++yzCA0N/UMRpN905HXevXs3PvroI2RkZHRBQsfFKyOCPPfcc5DJZJd8tPeXyO99++23+Pnnn/HOO+90bmgHZMvX+X/p9XpMnToVgwcPxt/+9rcrD07UhV5//XWsXr0aGzZsgLu7u+g4TqOqqgr33HMPli9fjoCAANFx7BqvjAjy5z//Gffee+8lj+nTpw+Cg4P/MDGqqakJFRUVFx1++fnnn5GVlQVfX99Wz992220YP348duzYcQXJHYstX+cWVVVVmDx5Mry9vbFhwwa4ubldaWynERAQAIVCgeLi4lbPFxcXX/R1DQ4Otup46tjr3GLJkiV4/fXX8dNPP2H48OG2jOnwrH2ds7KykJubi2nTplmeM5lMAMxXXk+fPo2+ffvaNrSjED1phS6tZWLlgQMHLM9t2bLlkhMri4qKpKNHj7Z6AJD+9a9/SdnZ2V0V3aF05HWWJEnS6XTS2LFjpQkTJkg1NTVdEdXhjBkzRvrTn/5k+W+j0SiFhYVdcgLrTTfd1Oq5+Ph4TmC9DGtfZ0mSpDfeeEPSaDRSampqV0R0Cta8znV1dX/4XTx9+nTpuuuuk44ePSo1NDR0ZXS7xjLiACZPnizFxsZK+/fvl3bv3i31799fmj17tuX/nz9/Xho4cKC0f//+i/4Z4N00l2Xt66zT6aS4uDhp2LBhUmZmplRUVGR5NDU1ifo27M7q1asltVotrVy5Ujpx4oT0wAMPSL6+vpJWq5UkSZLuuece6bnnnrMcv2fPHkmpVEpLliyRTp48KS1atEhyc3OTjh49KupbcAjWvs6vv/66pFKppHXr1rX62a2qqhL1LTgEa1/n3+PdNG1jGXEA5eXl0uzZs6Vu3bpJGo1Gmj9/fqtfGDk5ORIAafv27Rf9M1hGLs/a13n79u0SgDYfOTk5Yr4JO/Wf//xHioiIkFQqlTRmzBhp3759lv83YcIEad68ea2OX7t2rTRgwABJpVJJQ4YMkTZt2tTFiR2TNa9zr1692vzZXbRoUdcHdzDW/jz/L5aRtskkSZK6emiIiIiIqAXvpiEiIiKhWEaIiIhIKJYRIiIiEoplhIiIiIRiGSEiIiKhWEaIiIhIKJYRIiIiEoplhIiIiIRiGSEiIiKhWEaIiIhIKJYRIiIiEoplhIiIiIT6f9gfoxXs1GZJAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,sol,'-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br>\n",
    "<br>\n",
    "<br>\n",
    "Now we can proceed to compare the analytical solution with the numerical solution. <br>\n",
    "Remember, to execute the next cell you need to run first the sample utility."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "../postProcessing/sampleDict/20/s2_U.xy not found.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[12], line 6\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpyplot\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;28;01mas\u001b[39;00m\u001b[38;5;250m \u001b[39m\u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[1;32m      4\u001b[0m fig \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39mfigure(figsize\u001b[38;5;241m=\u001b[39m(\u001b[38;5;241m12\u001b[39m,\u001b[38;5;241m8\u001b[39m))\n\u001b[0;32m----> 6\u001b[0m data\u001b[38;5;241m=\u001b[39m\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mloadtxt\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m../postProcessing/sampleDict/20/s2_U.xy\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mskiprows\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m      8\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(data[:,\u001b[38;5;241m0\u001b[39m],data[:,\u001b[38;5;241m1\u001b[39m],\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mo\u001b[39m\u001b[38;5;124m'\u001b[39m,label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124micoFoam\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     10\u001b[0m vmax \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mmax\u001b[39m(data[:,\u001b[38;5;241m1\u001b[39m])\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/numpy/lib/npyio.py:1356\u001b[0m, in \u001b[0;36mloadtxt\u001b[0;34m(fname, dtype, comments, delimiter, converters, skiprows, usecols, unpack, ndmin, encoding, max_rows, quotechar, like)\u001b[0m\n\u001b[1;32m   1353\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(delimiter, \u001b[38;5;28mbytes\u001b[39m):\n\u001b[1;32m   1354\u001b[0m     delimiter \u001b[38;5;241m=\u001b[39m delimiter\u001b[38;5;241m.\u001b[39mdecode(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlatin1\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[0;32m-> 1356\u001b[0m arr \u001b[38;5;241m=\u001b[39m \u001b[43m_read\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcomment\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcomment\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdelimiter\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdelimiter\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1357\u001b[0m \u001b[43m            \u001b[49m\u001b[43mconverters\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mconverters\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mskiplines\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mskiprows\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43musecols\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43musecols\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1358\u001b[0m \u001b[43m            \u001b[49m\u001b[43munpack\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43munpack\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mndmin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mndmin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1359\u001b[0m \u001b[43m            \u001b[49m\u001b[43mmax_rows\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmax_rows\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mquote\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mquotechar\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1361\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m arr\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/numpy/lib/npyio.py:975\u001b[0m, in \u001b[0;36m_read\u001b[0;34m(fname, delimiter, comment, quote, imaginary_unit, usecols, skiplines, max_rows, converters, ndmin, unpack, dtype, encoding)\u001b[0m\n\u001b[1;32m    973\u001b[0m     fname \u001b[38;5;241m=\u001b[39m os\u001b[38;5;241m.\u001b[39mfspath(fname)\n\u001b[1;32m    974\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(fname, \u001b[38;5;28mstr\u001b[39m):\n\u001b[0;32m--> 975\u001b[0m     fh \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlib\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_datasource\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfname\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mrt\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    976\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m encoding \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    977\u001b[0m         encoding \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(fh, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mencoding\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlatin1\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/numpy/lib/_datasource.py:193\u001b[0m, in \u001b[0;36mopen\u001b[0;34m(path, mode, destpath, encoding, newline)\u001b[0m\n\u001b[1;32m    156\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    157\u001b[0m \u001b[38;5;124;03mOpen `path` with `mode` and return the file object.\u001b[39;00m\n\u001b[1;32m    158\u001b[0m \n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    189\u001b[0m \n\u001b[1;32m    190\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m    192\u001b[0m ds \u001b[38;5;241m=\u001b[39m DataSource(destpath)\n\u001b[0;32m--> 193\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mds\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpath\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnewline\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mnewline\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/numpy/lib/_datasource.py:533\u001b[0m, in \u001b[0;36mDataSource.open\u001b[0;34m(self, path, mode, encoding, newline)\u001b[0m\n\u001b[1;32m    530\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m _file_openers[ext](found, mode\u001b[38;5;241m=\u001b[39mmode,\n\u001b[1;32m    531\u001b[0m                               encoding\u001b[38;5;241m=\u001b[39mencoding, newline\u001b[38;5;241m=\u001b[39mnewline)\n\u001b[1;32m    532\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 533\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mFileNotFoundError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mpath\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m not found.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: ../postProcessing/sampleDict/20/s2_U.xy not found."
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 1200x800 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure(figsize=(12,8))\n",
    "\n",
    "data=np.loadtxt('../postProcessing/sampleDict/108/s2_U.xy', skiprows=0)\n",
    "\n",
    "plt.plot(data[:,0],data[:,1],'o',label='icoFoam')\n",
    "\n",
    "vmax = max(data[:,1])\n",
    "rmax = 0.5\n",
    "\n",
    "x = np.linspace(-1*rmax, rmax, 100)\n",
    "sol = vmax*(1 - x**2/rmax**2)\n",
    "\n",
    "plt.plot(x,sol,'-',label='Analytical solution')\n",
    "\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.xlabel('Radius',fontsize=15)\n",
    "plt.ylabel('$V_{axial}$',fontsize=15)\n",
    "\n",
    "\n",
    "\n",
    "#To save a figure\n",
    "#plt.savefig('test.png', format='png', dpi=300)\n",
    "#plt.savefig('test.pdf', format='pdf', dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For completeness, I will give you a few instructtions of how to work with arrays in numpy.<br>\n",
    "\n",
    "To create an array (in this case a vector):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.array([1, 2, 3, 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create a multi-dimensional array and save the values in the variable x:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To know the type of a variable you use type():"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "type(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which states that the variable x is a numpy ndarray. <br>\n",
    "To know the dimensions of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To print one value of the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x[0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To print one column of the array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and to print a row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x[0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the mean of the previous row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mean(x[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
